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The present study identifies a process for performing computational fluid dynamic 
calculations of the flow over full three-dimensional (3D) representations of complex ice 
shapes deposited on aircraft surfaces. Rime and glaze icing geometries formed on a 
NACA23012 airfoil were obtained during testing in the NASA Glenn Research Center’s 
Icing Research Tunnel (IRT). The ice shape geometries were scanned as a cloud of data 
points using a 3D laser scanner. The data point clouds were meshed using Geomagic 
software to create highly accurate models of the ice surface. The surface data was imported 
into Pointwise grid generation software to create the CFD surface and volume grids. It was 
determined that generating grids in Pointwise for complex 3D icing geometries was possible 
using various techniques that depended on the ice shape. Computations of the flow fields 
over these ice shapes were performed using the NASA National Combustion Code (NCC). 
Results for a rime ice shape for angle of attack conditions ranging from 0 to 10 degrees and 
for freestream Mach numbers of 0.10 and 0.18 are presented. For validation of the 
computational results, comparisons were made to test results from rapid-prototype models 
of the selected ice accretion shapes, obtained from a separate study in a subsonic wind tunnel 
at the University of Illinois at Urbana-Champaign. The computational and experimental 
results were compared for values of pressure coefficient and lift. Initial results show fairly 
good agreement for rime ice accretion simulations across the range of conditions examined. 
The glaze ice results are promising but require some further examination. 
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;; = distance above the local surface 

= non-dimensional wall distance for a wall-bounded flow 
- wall shear stress 
p = density of air 


I. Introduction 

Icing continues to be a major safety concern for the aircraft industry. Previous studies have shown that 
aerodynamic performance significantly degrades as ice accretes onto the wings. Depending on the ice shape 
geometry, there can be a substantial change in performance characteristics such as lift, drag, and pitching moment. 
Icing decreases the stall angle of attack, affecting the maximum coefficient of lift, while drastically increasing the 
drag of the wing which can lead to catastrophic results if the pilot is unaware of such effects. Therefore, icing 
research has been conducted extensively throughout the past several decades. 

Understanding the impact on the aerodynamic performance of aircraft wings resulting from changes to the 
surface geometry during in-flight ice deposition continues to be a topic of interest from the perspective of 
experimental and computational investigations. Although numerous experimental studies have been performed to 
obtain aerodynamic data for full 3D ice shape geometries^’^’^, computational studies to date have been limited to 2D 
cross-sections^’^, 3D extrusions of two dimensional cross-sections^’^, and smooth 3D ice shapes generated by ice 
accretion codes. ^ The ability to calculate the flow over a more accurate representation of an actual ice accretion 
surface has the potential to provide more information for determination of how ice shapes alter wing aerodynamics 
and how ice shape roughness impacts convective heat transfer rates during the accretion process. These 
improvements in turn could lead to more accurate ice accretion simulation tools and a more robust analysis of the 
impact of an icing encounter. 

Recent improvements in laser scanning technology and the associated post -processing software have enabled the 
recording of detailed ice shape surface geometry information. That digital representation of a complete ice shape, 
including surface roughness characteristics, can potentially be used to create rapid-prototype models of ice shapes 
and may also be used in grid generation software to develop realistic surface grids for use in computational fluid 
dynamic software. The goal of this study is to investigate the use of such data for the latter purpose in order to 
develop a methodology for determination of the aerodynamic impact of realistic ice shape geometries. The results of 
these initial studies indicate that the geometry data can be captured using commercial grid generation tools and that 
realistic aerodynamic results can be obtained. As this is an initial investigation into the use of the scanner data, a 
comprehensive examination of the details of the computational results was not performed but rather such an effort 
should form the basis for future studies. 


II. Methods 

In order to computationally analyze the 3D ice shape geometries, detailed surface grids representing the ice 
shape geometry need to be generated. In a separate study conducted in the Icing Research Tunnel (IRT) at NASA 
Glenn Research Center, a 3D laser scanner was used to collect a cloud of data points that correspond to surface 
locations for highly detailed 3D geometries of rime and glaze ice shapes. Lee, et al.^^ provides a description of the 
process used to scan the data. Photographs of the actual ice shapes used in this study are shown in Figure 1. 

Once the cloud of data points has been acquired, the Geomagic software package was used to generate 
representations of the surface suitable for use in grid generation software. Using Geomagic, the data points were 
connected in a grid-like manner to create polygonal cells representing the ice surface that captured the overall 
structure and the small-scale roughness features of the ice accretion. Due to the line-of-sight nature of the 3D 
scanning method, small patches with no data points were often found around highly concave geometric features. 
Interpolation is needed at these locations and was accomplished using features found within Geomagic. The end 
result of this post-processing effort is a water-tight surface representation of the scanned ice shape. Results of the 
scanning and post-processing with Geomagic are shown in Figure 2. 

Taking the surface models generated from Geomagic, CFD surface and volume grids were generated around the 
ice accreted airfoil models. Pointwise^^, a CFD grid generation software package, was used to generate the CFD 
grids for the icing geometry models. The output of the Geomagic software can be imported into Pointwise as either a 
database surface or as a surface grid. If it is imported as a database, the surface grid is constructed from that database 
information. If the geometry is imported as a grid, then the surface grid generation step can be skipped and the 
volume grid can be generated. This is only recommended if the geometry is very complex and grid domains cannot 
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a) Rime Ice Accretion b) Glaze Ice Accretion 


Figure 1. Rime and Glaze Ice Accretions on a NACA 23012 Airfoil Model 

be created around the geometry from the database format, such as the glaze ice accretion shown below. (Domains 
are the terminology used in the Pointwise software to identify surfaces and other boundaries of volumes created by 
the user in the process of making the overall grid structure.) Importing the geometry directly as a grid prevents the 
user from changing the number of cells on the grid, making refinement or smoothing impossible. During this study, 
the rime ice shapes were imported as databases while the glaze ice shapes were imported as grids. Importing the grid 
for the latter case was done because it was found that importing the more complex glaze shape as a database led 
problems in accurately projecting the domains onto the surface, which in turn made creation of a block 
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(a subset of the overall volume grid) problematic. The surface grid representation of the roughness for the rime ice 
shape could be changed within the Pointwise software but similar modifications for the glaze ice shapes had to be 
modified using Geomagic and then re-imported into Pointwise. A surface grid is shown in Figure 3 for the rime ice 
accretion and a surface and volume grid is shown for the glaze ice accretion in Figure 4. 



Figure 3. Surface Mesh for Rime Ice Accretion. 



Figure 4. Surface and Volume Mesh for a Glaze Ice Accretion. 

When creating any surface or volume grid, the user can choose to build the grid as a structured or unstructured 
grid. Structured grids enable faster computation times because each node is tabulated in a way such that every 
adjacent node is known. This tabulation allows the CFD solver to quickly solve the flow field within each cell 
because the solver can easily access the information in adjacent cells. Structured grids, however, are more difficult 
to manipulate to ftilly capture complex geometries such as ice shapes. The generation of a structured grid on the 
rime ice shape was attempted but in order to accurately capture the geometry, the overall grid had more than 200 
million cells which was considered excessive. Therefore, for this study, an unstructured grid was used for the ice 


4 

American Institute of Aeronautics and Astronautics 





accretion surfaces. An unstructured grid is good for capturing a lot of the detail in a complex geometry but requires 
more computation time per node, as the node locations are not stored with respect to adjacent nodes. Figures 3 and 4 
provide examples of the unstructured grids used for this study. 

Grid refinement is an important part of the grid generation process as it directly affects the accuracy of the flow 
solution. Refinement and concentrating cells in regions of aerodynamic importance is absolutely necessary to hilly 
capture flow features such as boundary layer growth, vortex shedding, trailing edge wake, and flow separation. This 
may require refinement of the grid aher initial calculations have been performed to understand where refinements 
are needed. Adaptive mesh technology could also be used for such refinement efforts but that was beyond the scope 
of this study. The amount of refinement needed to hilly capture the boundary layer for a viscous flow relies on the y^ 
value as well as hee stream flow parameters. Assuming the Reynolds number, free stream velocity, and density are 
known, a y^ value can be chosen which allows the turbulence model to fully capture the boundary layer. This in turn 
provides an estimate to the grid spacing normal to a wall surface. The following procedure gives an initial estimate 
to the wall distance for the first cell above the wing surface. First the skin hiction is estimated using the expression 
horn Schlichting shown in Equation (1). 


Cf = [21ogio(^?es) - 0.65]-2-3 for Re^ < 10^ (1) 

Then the wall shear stress is computed using Equation (2). 

T„ = Cf*\pUj (2) 

Using the wall shear stress, the hiction velocity can be calculated using Equation (3). 


u 


* 



Finally, the wall distance of the cell nearest to the surface can be calculated using Equation (4). 


y = ^ 


(3) 


(4) 


The y value can then be used as the initial cell size normal to the surface. 

In addition to grid spacing at the surface, grid refinement is also needed to make sure that cell skewness is 
minimized. Highly skewed cells can affect the convergence of the grid and in a worst case scenario they can prevent 
convergence entirely. Using an 0-grid or C-grid for airfoils and wings generally prevents skewed cells, but for these 
complex geometries skewness needs to be monitored and evaluated near the ice shape where very sharp angles are 
present. Use of tools available in grid generation sohware to monitor such grid quality metrics is an important part 
of the grid development process. 

If the volume grid is generated as an unstructured grid, Pointwise offers a unique tool called T-Rex gridding that 
generates layers of tetrahedral cells which better capture the flow in the boundary layer region. The number of 
layers, the growth rate of these layers, and the initial cell size in the first layer can all be varied to generate different 
size tetrahedral layers for different flow conditions. An example of the T-Rex gridding method as applied to the rime 
ice configuration is shown in Figure 5. 



a) Without T-Rex b) With T-Rex 


Figure 5. Example of the T-Rex Grid Generation Applied to Rime Ice Configuration. 
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A second major goal of this study was to investigate the feasibility of generating CFD grids that would converge 
to a valid flow solution in a CFD solver. Initially, clean wing simulations were performed for validation of 
coefficient of pressure, coefficient of lift and coefficient of drag aerodynamic performance properties using the CFD 
software selected for this study, the National Combustion Code^^'^^ (NCC). Rime and glaze ice accreted wings were 
then run in NCC for simulation validation with wind tunnel data. 

NCC has been and continues to be developed at the NASA Glenn Research Center and is used to perform the 
steady and unsteady Reynolds-averaged Navier-Stokes (RANS, URANS) simulations for this numerical study. The 
unstructured formulation of the main solver enables the mixing of various types of 2D and 3D elements for more 
efficiently modeling the problems of interest. Currently, it accepts quadrilateral and triangular meshes for two- 
dimensional geometries and tetrahedral, prismatic and hexahedral meshes for three-dimensional geometries. The 
code also takes advantage of the modem parallel algorithm so that it can be executed on a wide range of computer 
platforms to date. The parallelization is done via automatic domain decomposition that subdivides the entire domain 
to the number of available processors (or workstations) on the network. The underlying message passing library can 
be either PVM (Parallel Virtual Machine) or MPI (Message Passing Interface). 

The main solver, Corsair-CCD, is a finite-volume code with an explicit, four-stage Runge-Kutta integration 
algorithm. A time preconditioning technique is used to enhance the numerical efficiency for low speed flow 
computation. Local time stepping and residual smoothing are used to accelerate the convergence. Two-equation (k- 
8) turbulence model is used in the solver with the option of either generalized wall ftinction or direct integration as 
the turbulent wall boundary conditions. Furthermore, in addition to the standard k-s model, a higher order non-linear 
treatment of the turbulence stress term is included for better prediction of the turbulent flows with possible rotation 
and swirling^^, e.g. iced airfoil models and rotor blades in static/dynamic stall or in local ice-induced separated flow. 

In recent years, the NCC software team has ftirther developed a turbulence modeling approach^ called 
partially resolved numerical simulation (PRNS), that is capable of extracting some of the large-scale unsteady 
features of turbulent flows at reasonable computational costs. This approach affords an intermediate resolution of 
turbulence scales relative to those of RANS and LES and has the characteristics of the very large-eddy simulation 
(VLES). The very large scales of turbulence are directly calculated, and the effects of the unresolved scales are 
accounted for by an eddy viscosity model that is evolved from state-of-the-art models used in the RANS approach. 
NCC was developed and is used primarily for the development of combustion technology but its advancements^^ 
made in the areas of turbulence modeling, numerical simulation and computing platforms also greatly facilitates the 
application to 3D iced aerodynamics CFD simulations. This numerical platform has demonstrated its improved 
capability and flexibility with many benchmark validation cases^\ modeling complex internal reacting flow and 
external unsteady flow separation in static stall. These features are the reason for its use in this study. 

NCC requires specification of the following boundary conditions; inlet, exit, and side wall. In this study, periodic 
boundary conditions were used for the side walls. Periodic boundary conditions allow for any span -wise flow that is 
at one side to be transferred over to the other side which simulates an “infinite” wingspan where there is no wall 
effect to perturb the flow. This is not quite the situation for an airfoil with ice on the leading edge as the ice shape 
geometry on one side of the airfoil is not exactly the same as the geometry on the opposite side. Future studies could 
be conducted to examine the impact of this boundary condition selection for the analysis. For the inlet and exit 
boundary conditions, free stream conditions such as velocity, pressure, density, and temperature are required. 
Different values are used to change the effective angle of attack or altitude or Mach number that the geometry will 
experience. Figure 6 shows all of the prescribed boundary conditions for a C-grid case. 


Velocity 

Inlet 


Pressure Outlet 



Pressure 

Outlet 


Outlet 


Periodic 



Outlet 



Front View 


Isometric View 


Figure 6. Prescribed Boundary Conditions Used for C-grid in NCC. 
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The output from NCC only provides values for local gauge pressure and no other aerodynamic parameter 
because NCC was designed for analyzing combustors. Therefore, it was necessary to develop an external method of 
calculating coefficient of pressure, lift, and drag using the local gauge pressure. First the results were imported into 
Tecplot 360^^, a data viewing software package. There, a macro split up the 3D results into 2D cross-sections and 
calculated the coefficient of pressure using Equation (5). 

= 1 , 2 ( 5 ) 

P V2P«U« ^ ’ 

Using the coefficient of pressure values along the surface, the area underneath the coefficient of pressure was 
computed to calculate the normal and tangential force that would be used to calculate the coefficient of lift and drag. 
A trapezoidal area was assumed as the integration method. This computation was written in Python and performed 
as a post processing task after the NCC calculation. At the time of this writing, the drag computations for that script 
were not complete so only lift value comparisons will be presented in the results section. 

III. Results 

Computational results needed to be validated by comparison with experimental results. The ice shape geometry 
data generated from the scanning process and used to generate the CFD grids were also used to create rapid- 
prototype models which were tested in a subsonic wind tunnel at the University of Illinois at Urbana-Champaign^^. 
Pressure measurements as well as lift, drag, and moment measurements were recorded and calculated for an angle of 
attack sweep at two different test conditions. These test conditions are speeds of Mach 0.10 and 0.18, with a 
respective Reynolds numbers of 1.0x10^ and 1.8x10^. Static pressure was slightly below atmospheric, at 
approximately 98,600 Pa with a static temperature of 294 K. A sweep of angle of attack included -9 to 18 degrees 
during the test program. The same test conditions for the wind tunnel were used when running the CFD simulations 
using NCC, with the exception that we limited ourselves to 0 to 10 degrees incidence angle. 

As mentioned, NCC can be run in steady-state and unsteady modes. The calculations for this study were 
performed in an unsteady mode due to the nature of the flow near the surface. While the main flow on the scale of 
the airfoil model chord would be steady- state, it was anticipated that the flow over the ice shapes, especially around 
roughness elements, would be unsteady. As such, the pressure coefficient and lift results presented herein reflect the 
instantaneous nature of the flow. This is especially evident in the results for the glaze ice shape described below. 
Time averaging of the results was not performed during this study but could provide greater agreement with 
experiment and could form an element of future studies. 








Computational Results 
Experimental Results 
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Figure 7. Pressure Coefficient Results for the Rime Ice Geometry, Mid-Span Location. 
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Coefficient of pressure at a zero angle of attack for the rime ice accretion case is shown in Figure 7. This 
coefficient of pressure is taken at a half-span location of the 3D geometry. Preliminary results show fair 
correspondence between the experimental data and computational results but some anomalies are present. The 
pressure coefficient curves have a difference between experimental and computational results, specifically on the 
lower surface of the airfoil. Upper surface pressure correlates very well with a small discrepancy at the trailing edge. 
Integrating the coefficient of pressure across the upper and lower surface yields a coefficient of lift of 0.1366, 
resulting in a 15.4% difference with experimental data. 

Pressure contour plots are shown in Figure 8. These plots show the results for a clean (a) and rime ice (b) 
configuration. The contours show that the rime ice accretion does not substantially disturb the flow over the airfoil 
model. More detailed evaluations of the flow near the surface will be required to reveal the changes to the air flow 
resulting from the rime ice shape and its associated roughness features. Such an evaluation is made possible as a 
result of being able to capture the details of the rime ice surface geometry. That analysis is beyond the scope of the 
current study. 



a) clean airfoil model b) airfoil model with rime ice accretion 


Figure 8. Gauge Pressure Contour for NACA 23012 Airfoil Model, AOA = 0°, Re = 1.0x10^, M = 0.1 

Similarly, the coefficient of pressure at zero angle of attack is shown for the glaze ice accretion case. This 
coefficient of pressure is taken at a half-span location as well. Referring to Figure 9, the pressure coefficient varies 
largely at the leading edge of the airfoil due to the complex and irregular geometry. These fluctuations may also be a 
result of the unsteady nature of the flow around the horn region. Further examination, using time averaging of 
instantaneous results, may be needed to determine if that is the reason for the discrepancy in this region. 


Computational Results 
Experimental Results 
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Figure 9. Pressure Coefficieut Results for the Glaze Ice Geometry, Mid-Spau Locatiou. 
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Downstream of the leading edge, it is seen that both computational upper and lower pressure curves are 
underestimated in comparison to the experimental data. Further examination of this case is required before 
meaninghil results can be discussed. Some features of the computation that should be mentioned are the flow 
separation that occurs around the horns. Previous studies^^ have shown that higher grid resolution in the shear layer 
dividing the separated flow from the main free stream flow results in more realistic results for geometries such as 
this one. It is anticipated that more work along those lines is needed for this computation. For the purposes of this 
study it was our intention to determine whether a complete representation of a glaze ice shape could be prepared in 
commercial grid generation software and if a CFD solution could be obtained. 

The glaze ice results indicate that further examination of that geometry is required to understand the differences 
between computation and experiment. The pressure contour results shown in Figure 10 show a significant change 
from the pressure contours shown in Figure 8a for a clean NACA 23012 airfoil model. These initial results suggest 
that the ice shape geometry plays a major role in the outcome and that the surface roughness may be less influential. 
However, as for the rime ice results, further analysis will be necessary to confirm such observations. 



Figure 10. Gauge Pressure Contour for NACA 23012 Airfoil Model with Glaze Ice Accretion, 

AOA = 0°, Re = 1.0x10*, M = 0.1 

Since the rime ice results at zero degrees AOA reasonably matched the experiments, computations at several 
other angle of attack conditions and at a second Reynolds number, Mach number combination were performed. The 
results are summarized in the lift versus angle of attack plots shown in Figure 11. 

The average percent difference between computation and experiment for the low Re, M case is 15.2% with a 
minimum difference of 6.8% and a maximum difference of 21.3%. For the higher Re, M condition, the average 
difference is 11.5% with a minimum difference of 8.1% and a maximum difference of 19.4%. The absolute 
difference between computation and experiment increases with increasing angle of attack for both the low Re, M 
condition as well as the high Re, M condition. 

These results indicate two trends. As the angle of attack increases, the difference between the computation and 
experiment increases. Additionally, as the Reynolds number and Mach number increase the differences between the 
computation and experiment decrease. With the angle of attack increase the influence of the ice shape on the 
pressure profile becomes more important. Thus, accurate modeling of the unsteady and locally separated flow 
becomes more important. The results suggest that further study is required to understand the influence of ice shape 
measurement resolution and associated surface grid spacing to better understand the requirements for accurate 
simulation. This also suggests that greater attention to the drag results should be considered in future work as 
detailed examination of the surface flow features should provide clues on how best to simulate these shapes. 
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a) Re = 1.0x10®, M = 0.1 


b) Re = 1.8x10®, M = 0.18 


Figure 11. Comparison of Computational and Experimental Lift Coefficient for Rime Ice Shape on 

NACA 23012 Airfoil Model. 

The closer agreement between computation and experiment at the higher Reynolds number and Mach number 
also indicate that accurate simulation of the viscous region plays a role in the overall computation. When the viscous 
influence on the results decreases as in the case shown here there becomes better agreement and the influence of the 
ice shape roughness details seems to have a lesser influence. These are, of course, just a single set of computations 
and hirther examination of the results as well as more comparisons to other ice shape results will be required to 
verify such trends. 


IV. Conclusion 

This study shows that it is possible to generate CFD grids representing highly complex 3D icing geometries and 
to perform viscous computational analysis on these CFD grids. Preliminary results show some correlation between 
computational results and experimental data but a general trend showed that computational pressure coefficient 
curves deviate somewhat from the experimental data. This discrepancy results in lift coefficient values that are 
larger than the values obtained from experiment. Further work needs to be done to determine what modifications to 
the analysis process will lead to higher correlation. Recommendations include running a parametric study on the 
refinement of the CFD grids, evaluate the turbulence models impact on the computations, develop a complete and 
robust method of calculating the 3D lift and drag forces from the computational surface pressures, and complete the 
angle of attack sweep for both the rime and glaze icing geometries to get a frill comparison with the experimental 
data. 

This study indicates that existing grid generation and CFD tools have the potential to provide research personnel 
with the ability to calculate the flow over a more accurate representation of an actual ice accretion surface. Such a 
capability has the potential to provide more information for determination of how ice shapes alter wing 
aerodynamics and how ice shape roughness impacts convective heat transfer rates during the accretion process. 
These improvements in turn could lead to more accurate ice accretion simulation tools and a more robust analysis of 
the impact of an icing encounter. 
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